#script for mixed model analysis of Ek and irradiance

Load the various libraries

library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(car)
library(MuMIn)
library (dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(emmeans)
library(DHARMa)
## This is DHARMa 0.4.5. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)
library(patchwork)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ tibble  3.1.8     ✔ purrr   0.3.5
## ✔ tidyr   1.2.1     ✔ stringr 1.4.1
## ✔ readr   2.1.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()   masks Matrix::expand()
## ✖ rstatix::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ tidyr::pack()     masks Matrix::pack()
## ✖ dplyr::recode()   masks car::recode()
## ✖ purrr::some()     masks car::some()
## ✖ tidyr::unpack()   masks Matrix::unpack()
library(sjPlot)
library(sjmisc)
## 
## Attaching package: 'sjmisc'
## 
## The following object is masked from 'package:purrr':
## 
##     is_empty
## 
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## 
## The following object is masked from 'package:tibble':
## 
##     add_case
library(mmtable2)
## 
## Attaching package: 'mmtable2'
## 
## The following object is masked from 'package:tidyr':
## 
##     table1
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load the data

ek_irrad_data <- read.csv("input_data/mean_supersaturation_by_period.csv")

Assign run as a factor

ek_irrad_data$Run <- as.factor(ek_irrad_data$Run)

Assign temperature as a factor

ek_irrad_data$Temperature <- as.factor(ek_irrad_data$Temp...C.)

Assigns treatment as characters from integers then to factors

ek_irrad_data$Treatment <- as.factor(as.character(ek_irrad_data$Treatment))

Toggle between the species for output. Use Day 9 for final analysis recent change: removing the very odd ek.1 value of 559.4 in hypnea dataset

hypnea <- subset(ek_irrad_data, Species == "hm")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

ulva <- subset(ek_irrad_data, Species == "ul")
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

Add growth rate from other dataset to this one and subset by species

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_102222.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)

Subset the growth rate data by Species

gr_ulva <- subset(growth_rate, Species == "Ul")
gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017)

Calculate the growth rate from the data

ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

#_______________________ULVA____________________________

Run model without interaction between the treatments and temperature - irradiance_over_ek is in minutes. Take RLC.Order out of the random effects because it is causing problems of singularity. R2 is same with or without (+ (1 | RLC.Order))

ek_irrad_model_ulva <- lmer(formula = supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = ulva)

Get an r2 for the model

r.squaredGLMM(ek_irrad_model_ulva)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.3642221 0.8897065

Make histograms of data and residuals (first test for mixed model use) and residual plots of the data for ulva

hist(ulva$supersat_total, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)

hist(resid(ek_irrad_model_ulva))

plot(resid(ek_irrad_model_ulva) ~ fitted(ek_irrad_model_ulva))

qqnorm(resid(ek_irrad_model_ulva))
qqline(resid(ek_irrad_model_ulva))

ulva %>% ggplot(aes(supersat_total)) +
        geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Check the performance of the model and run a summary

performance::check_model(ek_irrad_model_ulva)

summary(ek_irrad_model_ulva)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: ulva
## 
## REML criterion at convergence: 2813
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1222 -0.4370  0.0455  0.4915  2.8439 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept)  655.3   25.60   
##  Run      (Intercept) 2719.5   52.15   
##  Residual              708.3   26.61   
## Number of obs: 288, groups:  Plant.ID, 144; Run, 8
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    477.5252    37.4668    5.0114  12.745 5.21e-05 ***
## Treatment1     -50.3702    44.2995    4.9966  -1.137   0.3071    
## Treatment2     -47.6295    44.2995    4.9966  -1.075   0.3315    
## Treatment2.5  -157.9977    64.3119    4.8337  -2.457   0.0592 .  
## Treatment3     -80.6781    44.2995    4.9966  -1.821   0.1283    
## Treatment4     -88.4096    44.2995    4.9966  -1.996   0.1025    
## Temperature27    0.1846     6.8448  115.6816   0.027   0.9785    
## Temperature30   -2.3436     6.8448  115.6816  -0.342   0.7327    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.836                                          
## Treatment2  -0.836  0.992                                   
## Treatmnt2.5 -0.576  0.487  0.487                            
## Treatment3  -0.836  0.992  0.992  0.487                     
## Treatment4  -0.836  0.992  0.992  0.487  0.992              
## Temperatr27 -0.091  0.000  0.000  0.000  0.000  0.000       
## Temperatr30 -0.091  0.000  0.000  0.000  0.000  0.000  0.500

Check for equal variance to determine which ANOVA to use

bartlett.test(supersat_total ~ Treatment, data = ulva)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  supersat_total by Treatment
## Bartlett's K-squared = 19.346, df = 5, p-value = 0.001656

Run Welch’s ANOVA if not equal variance. Code forces one independent variable at a time. Games Howell Test for pairwise comparison for the only independent variable that was found to have significant results.

welch_anova_treatment_ulva <- oneway.test(supersat_total ~ Treatment, data = ulva, var.equal = FALSE)
welch_anova_treatment_ulva
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_total and Treatment
## F = 74.182, num df = 5.00, denom df = 130.53, p-value < 2.2e-16
welch_anova_temp_ulva <- oneway.test(supersat_total ~ Temperature, data = ulva, var.equal = FALSE)
welch_anova_temp_ulva
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_total and Temperature
## F = 0.16463, num df = 2.00, denom df = 189.68, p-value = 0.8483
games_howell_test(ulva, supersat_total ~ Treatment, conf.level = 0.95, detailed = FALSE)
## # A tibble: 15 × 8
##    .y.            group1 group2 estimate conf.low conf.high    p.adj p.adj.sig…¹
##  * <chr>          <chr>  <chr>     <dbl>    <dbl>     <dbl>    <dbl> <chr>      
##  1 supersat_total 0      1        -44.9     -77.4   -12.4   2   e- 3 **         
##  2 supersat_total 0      2        -42.1     -74.8    -9.53  4   e- 3 **         
##  3 supersat_total 0      2.5     -158.     -183.   -133.    4.75e-10 ****       
##  4 supersat_total 0      3        -75.2    -108.    -42.2   4.24e- 8 ****       
##  5 supersat_total 0      4        -82.9    -115.    -50.4   1.34e- 9 ****       
##  6 supersat_total 1      2          2.74    -34.4    39.9   1   e+ 0 ns         
##  7 supersat_total 1      2.5     -113.     -144.    -82.2   0        ****       
##  8 supersat_total 1      3        -30.3     -67.8     7.23  1.85e- 1 ns         
##  9 supersat_total 1      4        -38.0     -75.1    -0.939 4.1 e- 2 *          
## 10 supersat_total 2      2.5     -116.     -147.    -84.8   0        ****       
## 11 supersat_total 2      3        -33.0     -70.7     4.59  1.19e- 1 ns         
## 12 supersat_total 2      4        -40.8     -78.0    -3.58  2.3 e- 2 *          
## 13 supersat_total 2.5    3         82.8      51.3   114.    5.32e-10 ****       
## 14 supersat_total 2.5    4         75.1      44.1   106.    7.93e- 9 ****       
## 15 supersat_total 3      4         -7.73    -45.3    29.8   9.91e- 1 ns         
## # … with abbreviated variable name ¹​p.adj.signif
plot(allEffects(ek_irrad_model_ulva))

Summarize the means for rETRmax

ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_total), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          477.
## 2 1          432.
## 3 2          435.
## 4 2.5        319.
## 5 3          402.
## 6 4          394.
ulva %>% group_by(Run) %>% summarise_at(vars(supersat_total), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      461.
## 2 2      465.
## 3 3      422.
## 4 3.5    352.
## 5 4      349.
## 6 6      452.
## 7 6.5    501.
## 8 9      319.

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_irrad_ek_graph <- ggplot(ulva, aes(x=supersat_total, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "A", subtitle = "Ulva lactuca", x = "Mean saturation time (min)", 
             y = "growth rate (%)") + stat_regline_equation(label.x = 450, label.y = 225) + 
        stat_cor(label.x = 450, label.y = 215) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
ulva_growth_irrad_ek_graph
## `geom_smooth()` using formula 'y ~ x'

Plot a histogram in ggplot2

ulva %>% ggplot(aes(supersat_total)) +
        geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Make a boxplot for publications of the results

ulva %>% ggplot(aes(treatment_graph, supersat_total)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Mean Daily Saturation Time (min)", title= "A", subtitle = "Ulva lactuca") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(100, 600) + stat_mean() + 
        geom_hline(yintercept=100, color = "orange", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

##Run model for relative Hsat

rel_hsat_model_ulva <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = ulva)
hist(ulva$supersat_rel, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)

hist(resid(rel_hsat_model_ulva))

plot(resid(rel_hsat_model_ulva) ~ fitted(rel_hsat_model_ulva))

qqnorm(resid(rel_hsat_model_ulva))
qqline(resid(rel_hsat_model_ulva))

ulva %>% ggplot(aes(supersat_total)) +
        geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Check the performance of the model

performance::check_model(rel_hsat_model_ulva)

r.squaredGLMM(rel_hsat_model_ulva)
##            R2m       R2c
## [1,] 0.4420469 0.8743361
summary(rel_hsat_model_ulva)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: ulva
## 
## REML criterion at convergence: -802
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2462 -0.4333  0.0519  0.4954  2.7114 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept) 0.001543 0.03928 
##  Run      (Intercept) 0.004649 0.06818 
##  Residual             0.001800 0.04242 
## Number of obs: 288, groups:  Plant.ID, 144; Run, 8
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     0.751352   0.049313   5.008835  15.236 2.18e-05 ***
## Treatment1     -0.083563   0.058292   4.988804  -1.434   0.2113    
## Treatment2     -0.079353   0.058292   4.988804  -1.361   0.2317    
## Treatment2.5   -0.260459   0.084333   4.760091  -3.088   0.0290 *  
## Treatment3     -0.131558   0.058292   4.988804  -2.257   0.0738 .  
## Treatment4     -0.143701   0.058292   4.988804  -2.465   0.0570 .  
## Temperature27   0.009256   0.010659 112.449439   0.868   0.3870    
## Temperature30   0.010096   0.010659 112.449439   0.947   0.3456    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.833                                          
## Treatment2  -0.833  0.989                                   
## Treatmnt2.5 -0.576  0.487  0.487                            
## Treatment3  -0.833  0.989  0.989  0.487                     
## Treatment4  -0.833  0.989  0.989  0.487  0.989              
## Temperatr27 -0.108  0.000  0.000  0.000  0.000  0.000       
## Temperatr30 -0.108  0.000  0.000  0.000  0.000  0.000  0.500

Check for equal variance

bartlett.test(supersat_rel ~ Treatment, data = ulva)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  supersat_rel by Treatment
## Bartlett's K-squared = 10.634, df = 5, p-value = 0.05915
#run Welch's ANOVA if not equal variance
welch_anova_treatment_ulva <- oneway.test(supersat_rel ~ Treatment, data = ulva, var.equal = FALSE)
welch_anova_treatment_ulva
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_rel and Treatment
## F = 84.949, num df = 5.00, denom df = 131.02, p-value < 2.2e-16
welch_anova_temp_ulva <- oneway.test(supersat_rel ~ Temperature, data = ulva, var.equal = FALSE)
welch_anova_temp_ulva
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_rel and Temperature
## F = 0.5808, num df = 2.00, denom df = 189.75, p-value = 0.5604
games_howell_test(ulva, supersat_rel ~ Treatment, conf.level = 0.95, detailed = TRUE)
## # A tibble: 15 × 14
##    .y.       group1 group2    n1    n2 estimate conf.…¹ conf.h…²      se stati…³
##  * <chr>     <chr>  <chr>  <int> <int>    <dbl>   <dbl>    <dbl>   <dbl>   <dbl>
##  1 supersat… 0      1         48    48 -0.0770  -0.124  -3.06e-2 0.0113    4.83 
##  2 supersat… 0      2         48    48 -0.0728  -0.119  -2.63e-2 0.0113    4.56 
##  3 supersat… 0      2.5       48    48 -0.260   -0.299  -2.22e-1 0.00943  19.5  
##  4 supersat… 0      3         48    48 -0.125   -0.173  -7.72e-2 0.0116    7.62 
##  5 supersat… 0      4         48    48 -0.137   -0.184  -9.01e-2 0.0114    8.49 
##  6 supersat… 1      2         48    48  0.00421 -0.0463  5.47e-2 0.0123    0.243
##  7 supersat… 1      2.5       48    48 -0.183   -0.227  -1.40e-1 0.0106   12.3  
##  8 supersat… 1      3         48    48 -0.0480  -0.0996  3.66e-3 0.0126    2.70 
##  9 supersat… 1      4         48    48 -0.0601  -0.111  -9.15e-3 0.0124    3.43 
## 10 supersat… 2      2.5       48    48 -0.188   -0.231  -1.44e-1 0.0106   12.5  
## 11 supersat… 2      3         48    48 -0.0522  -0.104  -4.81e-4 0.0126    2.94 
## 12 supersat… 2      4         48    48 -0.0643  -0.115  -1.33e-2 0.0124    3.67 
## 13 supersat… 2.5    3         48    48  0.135    0.0904  1.80e-1 0.0109    8.77 
## 14 supersat… 2.5    4         48    48  0.123    0.0790  1.68e-1 0.0107    8.12 
## 15 supersat… 3      4         48    48 -0.0121  -0.0643  4.00e-2 0.0127    0.677
## # … with 4 more variables: df <dbl>, p.adj <dbl>, p.adj.signif <chr>,
## #   method <chr>, and abbreviated variable names ¹​conf.low, ²​conf.high,
## #   ³​statistic
plot(allEffects(rel_hsat_model_ulva))

ulva %>% ggplot(aes(treatment_graph, supersat_rel)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Mean Rel. Hsat Time", title= "A", subtitle = "Ulva lactuca") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(0, 1) + stat_mean() + 
        geom_hline(yintercept=0.5, color = "orange", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

Summarize the means for relative Hsat

ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_rel), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0         0.758
## 2 1         0.681
## 3 2         0.685
## 4 2.5       0.497
## 5 3         0.633
## 6 4         0.621
ulva %>% group_by(Run) %>% summarise_at(vars(supersat_rel), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1     0.704
## 2 2     0.723
## 3 3     0.674
## 4 3.5   0.567
## 5 4     0.572
## 6 6     0.723
## 7 6.5   0.792
## 8 9     0.497
ulva %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      643.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 6      626.
## 7 6.5    633.
## 8 9      641.
ulva_growth_rel_hsat_graph <- ggplot(ulva, aes(x=supersat_rel, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "A", subtitle = "Ulva lactuca", x = "Mean Rel. Hsat", 
             y = "growth rate (%)") + stat_regline_equation(label.x = 0.68, label.y = 200) + 
        stat_cor(label.x = 0.68, label.y = 215) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
ulva_growth_rel_hsat_graph
## `geom_smooth()` using formula 'y ~ x'

#______________________HYPNEA___________________________

Run model for Hypnea

ek_irrad_model_hypnea <- lmer(formula = supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)

Get an r2 for the model

r.squaredGLMM(ek_irrad_model_hypnea)
##            R2m       R2c
## [1,] 0.1460739 0.7597093

Make histograms of data and residuals (first test for mixed model use) and residual plots of the data for Hypnea

hist(hypnea$supersat_total, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hist(resid(ek_irrad_model_hypnea))

plot(resid(ek_irrad_model_hypnea) ~ fitted(ek_irrad_model_hypnea))

qqnorm(resid(ek_irrad_model_hypnea))
qqline(resid(ek_irrad_model_hypnea))

Plot a histogram in ggplot

hypnea %>% ggplot(aes(supersat_total)) +
        geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Check the performance of the model

performance::check_model(ek_irrad_model_hypnea)

Run a summary of the model

summary(ek_irrad_model_hypnea)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: hypnea
## 
## REML criterion at convergence: 2904.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3115 -0.4540  0.0553  0.5465  2.1801 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept)  787.4   28.06   
##  Run      (Intercept) 2032.5   45.08   
##  Residual             1104.2   33.23   
## Number of obs: 287, groups:  Plant.ID, 144; Run, 8
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    396.535     32.811   5.043  12.086 6.47e-05 ***
## Treatment1     -49.539     38.782   5.021  -1.277 0.257345    
## Treatment2     -35.766     38.782   5.021  -0.922 0.398558    
## Treatment2.5   -37.009     55.925   4.729  -0.662 0.538991    
## Treatment3     -51.386     38.782   5.021  -1.325 0.242263    
## Treatment4     -69.747     38.789   5.025  -1.798 0.131789    
## Temperature27   29.394      7.910 131.267   3.716 0.000299 ***
## Temperature30   32.267      7.915 131.714   4.077 7.86e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.830                                          
## Treatment2  -0.830  0.985                                   
## Treatmnt2.5 -0.575  0.487  0.487                            
## Treatment3  -0.830  0.985  0.985  0.487                     
## Treatment4  -0.830  0.985  0.985  0.487  0.985              
## Temperatr27 -0.121  0.000  0.000  0.000  0.000  0.000       
## Temperatr30 -0.121  0.000  0.000  0.000  0.000  0.001  0.500

Check for equal variance

bartlett.test(supersat_total ~ Treatment, data = hypnea)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  supersat_total by Treatment
## Bartlett's K-squared = 50.954, df = 5, p-value = 8.84e-10

Run Welch’s ANOVA if not equal variance. Code forces one independent variable at a time. Games Howell Test for pairwise comparison for the only independent variable that was found to have significant results.

welch_anova_treatment <- oneway.test(supersat_total ~ Treatment, data = hypnea, var.equal = FALSE)
welch_anova_treatment
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_total and Treatment
## F = 7.4774, num df = 5.00, denom df = 127.19, p-value = 3.509e-06
welch_anova_temp <- oneway.test(supersat_total ~ Temperature, data = hypnea, var.equal = FALSE)
welch_anova_temp
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_total and Temperature
## F = 11.458, num df = 2.00, denom df = 184.95, p-value = 2.036e-05
games_howell_test(hypnea, supersat_total ~ Treatment, conf.level = 0.95, detailed = FALSE)
## # A tibble: 15 × 8
##    .y.            group1 group2 estimate conf.low conf.high      p.adj p.adj.s…¹
##  * <chr>          <chr>  <chr>     <dbl>    <dbl>     <dbl>      <dbl> <chr>    
##  1 supersat_total 0      1        -42.7     -81.9     -3.53 0.025      *        
##  2 supersat_total 0      2        -28.9     -63.5      5.60 0.153      ns       
##  3 supersat_total 0      2.5      -37.0     -60.0    -14.0  0.000153   ***      
##  4 supersat_total 0      3        -44.6     -79.3     -9.82 0.004      **       
##  5 supersat_total 0      4        -62.0     -94.5    -29.6  0.00000416 ****     
##  6 supersat_total 1      2         13.8     -30.4     57.9  0.944      ns       
##  7 supersat_total 1      2.5        5.71    -30.6     42.0  0.997      ns       
##  8 supersat_total 1      3         -1.85    -46.2     42.5  1          ns       
##  9 supersat_total 1      4        -19.3     -61.9     23.3  0.773      ns       
## 10 supersat_total 2      2.5       -8.07    -39.2     23.1  0.973      ns       
## 11 supersat_total 2      3        -15.6     -56.0     24.7  0.869      ns       
## 12 supersat_total 2      4        -33.1     -71.5      5.37 0.134      ns       
## 13 supersat_total 2.5    3         -7.55    -39.0     23.8  0.98       ns       
## 14 supersat_total 2.5    4        -25.0     -53.8      3.79 0.125      ns       
## 15 supersat_total 3      4        -17.5     -56.1     21.2  0.776      ns       
## # … with abbreviated variable name ¹​p.adj.signif
games_howell_test(hypnea, supersat_total ~ Temperature, conf.level = 0.95, detailed = FALSE)
## # A tibble: 3 × 8
##   .y.            group1 group2 estimate conf.low conf.high     p.adj p.adj.sig…¹
## * <chr>          <chr>  <chr>     <dbl>    <dbl>     <dbl>     <dbl> <chr>      
## 1 supersat_total 20     27        38.3      16.4      60.2 0.000161  ***        
## 2 supersat_total 20     30        43.3      21.0      65.6 0.0000254 ****       
## 3 supersat_total 27     30         4.96    -12.7      22.6 0.784     ns         
## # … with abbreviated variable name ¹​p.adj.signif

If ANOVA more appropriate (not this case) run from code above with Tukey’s HSD pairwise comparisons where appropriate

Plot results for basic output

plot(allEffects(ek_irrad_model_hypnea))

Plot a regression between the photosynthetic independent variables of interest and growth rate and a histogram of the data

hypnea_growth_irrad_ek_graph <- ggplot(hypnea, aes(x=supersat_total, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "B", subtitle = "Hypnea musciformis", x = "Mean saturation time (min)", y = "growth rate (%)") + 
        stat_regline_equation(label.x = 400, label.y = 150) + stat_cor(label.x = 400, label.y = 140) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_irrad_ek_graph
## `geom_smooth()` using formula 'y ~ x'

Make a boxplot for publications of the results

hypnea %>% ggplot(aes(treatment_graph, supersat_total)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Mean Daily Saturation Time (min)", title= "B", subtitle = "Hypnea musciformis") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(100, 600) + stat_mean() + 
        geom_hline(yintercept=100, color = "orange", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
        theme_bw() +
        theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

Summarize the means for rETRmax

hypnea %>% group_by(Treatment) %>% summarise_at(vars(supersat_total), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          417.
## 2 1          374.
## 3 2          388.
## 4 2.5        380.
## 5 3          373.
## 6 4          355.
hypnea %>% group_by(Run) %>% summarise_at(vars(supersat_total), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      407.
## 2 2      426.
## 3 3      360.
## 4 3.5    312.
## 5 4      320.
## 6 7      380.
## 7 8      406.
## 8 9      428.

##Run the model for Hsat

rel_hsat_model_hypnea  <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)

Make a histogram and residual plots of the data for ulva

hist(hypnea$supersat_rel, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hist(resid(rel_hsat_model_hypnea))

plot(resid(rel_hsat_model_hypnea) ~ fitted(rel_hsat_model_hypnea))

qqnorm(resid(rel_hsat_model_hypnea))
qqline(resid(rel_hsat_model_hypnea))

hypnea %>% ggplot(aes(supersat_total)) +
        geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Check the performance of the model

performance::check_model(rel_hsat_model_hypnea)

r.squaredGLMM(rel_hsat_model_hypnea)
##            R2m       R2c
## [1,] 0.2080092 0.7254956
summary(rel_hsat_model_hypnea)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: hypnea
## 
## REML criterion at convergence: -701.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4632 -0.4603  0.0425  0.5323  2.3552 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept) 0.001822 0.04269 
##  Run      (Intercept) 0.003371 0.05806 
##  Residual             0.002755 0.05249 
## Number of obs: 287, groups:  Plant.ID, 144; Run, 8
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     0.61738    0.04279   5.05849  14.429 2.65e-05 ***
## Treatment1     -0.07498    0.05056   5.02872  -1.483   0.1978    
## Treatment2     -0.05364    0.05056   5.02872  -1.061   0.3369    
## Treatment2.5   -0.10415    0.07243   4.61640  -1.438   0.2147    
## Treatment3     -0.07807    0.05056   5.02872  -1.544   0.1828    
## Treatment4     -0.10692    0.05057   5.03413  -2.114   0.0878 .  
## Temperature27   0.05385    0.01223 133.00960   4.402 2.18e-05 ***
## Temperature30   0.06321    0.01224 133.48716   5.164 8.60e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.823                                          
## Treatment2  -0.823  0.978                                   
## Treatmnt2.5 -0.575  0.486  0.486                            
## Treatment3  -0.823  0.978  0.978  0.486                     
## Treatment4  -0.823  0.977  0.977  0.486  0.977              
## Temperatr27 -0.143  0.000  0.000  0.000  0.000  0.000       
## Temperatr30 -0.143  0.000  0.000  0.000  0.000  0.001  0.500

Check for equal variance

bartlett.test(supersat_rel ~ Treatment, data = hypnea)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  supersat_rel by Treatment
## Bartlett's K-squared = 49.512, df = 5, p-value = 1.744e-09
#run Welch's ANOVA if not equal variance
welch_anova_treatment <- oneway.test(supersat_rel ~ Treatment, data = hypnea, var.equal = FALSE)
welch_anova_treatment
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_rel and Treatment
## F = 16.025, num df = 5.00, denom df = 127.18, p-value = 3.143e-12
welch_anova_temp <- oneway.test(supersat_rel ~ Temperature, data = hypnea, var.equal = FALSE)
welch_anova_temp
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  supersat_rel and Temperature
## F = 16.829, num df = 2.00, denom df = 185.09, p-value = 1.927e-07
games_howell_test(hypnea, supersat_rel ~ Treatment, conf.level = 0.95, detailed = TRUE)
## # A tibble: 15 × 14
##    .y.       group1 group2    n1    n2 estimate conf.…¹ conf.h…²      se stati…³
##  * <chr>     <chr>  <chr>  <int> <int>    <dbl>   <dbl>    <dbl>   <dbl>   <dbl>
##  1 supersat… 0      1         48    48 -0.0663  -0.126  -0.00664 0.0144    3.25 
##  2 supersat… 0      2         48    48 -0.0449  -0.0959  0.00601 0.0124    2.57 
##  3 supersat… 0      2.5       48    48 -0.104   -0.140  -0.0687  0.00858   8.58 
##  4 supersat… 0      3         48    48 -0.0694  -0.122  -0.0171  0.0127    3.87 
##  5 supersat… 0      4         48    47 -0.0973  -0.146  -0.0488  0.0118    5.85 
##  6 supersat… 1      2         48    48  0.0213  -0.0440  0.0866  0.0159    0.952
##  7 supersat… 1      2.5       48    48 -0.0379  -0.0926  0.0168  0.0131    2.04 
##  8 supersat… 1      3         48    48 -0.00309 -0.0694  0.0632  0.0161    0.136
##  9 supersat… 1      4         48    47 -0.0310  -0.0944  0.0324  0.0154    1.42 
## 10 supersat… 2      2.5       48    48 -0.0592  -0.104  -0.0143  0.0108    3.87 
## 11 supersat… 2      3         48    48 -0.0244  -0.0832  0.0343  0.0143    1.21 
## 12 supersat… 2      4         48    47 -0.0523  -0.108   0.00310 0.0135    2.75 
## 13 supersat… 2.5    3         48    48  0.0348  -0.0117  0.0812  0.0112    2.20 
## 14 supersat… 2.5    4         48    47  0.00689 -0.0351  0.0489  0.0101    0.482
## 15 supersat… 3      4         48    47 -0.0279  -0.0845  0.0287  0.0138    1.43 
## # … with 4 more variables: df <dbl>, p.adj <dbl>, p.adj.signif <chr>,
## #   method <chr>, and abbreviated variable names ¹​conf.low, ²​conf.high,
## #   ³​statistic
games_howell_test(hypnea, supersat_rel ~ Temperature, conf.level = 0.95, detailed = TRUE)
## # A tibble: 3 × 14
##   .y.    group1 group2    n1    n2 estim…¹ conf.…² conf.…³      se stati…⁴    df
## * <chr>  <chr>  <chr>  <int> <int>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>
## 1 super… 20     27        96    96  0.0674  0.0350  0.0999 0.00971   4.91   166.
## 2 super… 20     30        96    95  0.0785  0.0455  0.112  0.00989   5.62   171.
## 3 super… 27     30        96    95  0.0111 -0.0152  0.0374 0.00786   0.998  188.
## # … with 3 more variables: p.adj <dbl>, p.adj.signif <chr>, method <chr>, and
## #   abbreviated variable names ¹​estimate, ²​conf.low, ³​conf.high, ⁴​statistic
plot(allEffects(rel_hsat_model_hypnea))

hypnea %>% ggplot(aes(treatment_graph, supersat_rel)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Mean Rel. Hsat Time", title= "B", subtitle = "Hypnea musciformis") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(0, 1) + stat_mean() + 
        geom_hline(yintercept=0.5, color = "orange", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
        theme_bw() +
        theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

hypnea_growth_rel_hsat_graph <- ggplot(hypnea, aes(x=supersat_rel, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "B", subtitle = "Hypnea musciformis", x = "Mean Hsat Time (min)", y = "growth rate (%)") + 
        stat_regline_equation(label.x = 0.68, label.y = 100) + stat_cor(label.x = 0.68, label.y = 108) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_rel_hsat_graph
## `geom_smooth()` using formula 'y ~ x'